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Abstract 

The Hf-O system has been modeled by combining existing experimental data and first-principles calculations results through the 
CALPHAD approach. Special quasirandom structures of a and j3 hafnium were generated to calculate the mixing behavior of 
oxygen and vacancies. For the total energy of oxygen, vibrational, rotational and translational degrees of freedom were considered. 
The Hf-O system was combined with previously modeled Hf-Si and Si-O systems, and the ternary compound in the Hf-Si-O system, 
HfSiCU has been introduced to calculate the stability diagrams pertinent to the thin film processing. 
© 2008 Elsevier Ltd. All rights reserved. 
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1. Introduction 

The Hf-O system has been considered as one of the most 
important systems in various industrial fields, such as nu- 
clear materials and high temperature/pressure materials. 
Hafnium dioxide is known as the least volatile of all the 
oxides and its high melting point, extreme chemical inert- 
ness, and high thermal neutron capture cross section make 
it suitable for use as control rods or neutron shielding. All 
of these reasons can make Hf02 a promising refractory ma- 
terial for future nuclear applications [1]. 

Recently, the Hf-0 system has even attracted consider- 
able attention for semi-conductor materials. The current 
gate oxide material (Si02 in general) thickness in advanced 
complementary metal-oxide semiconductor (CMOS) inte- 
grated circuits has continuously decreased and has reached 
the current process limit [2]. One solution to further im- 
prove their performance is to use alternative materials with 
higher dielectric constants (k) , such as Zr02 and Hf02 [3, 4] . 
In this regard, thermodynamic stability calculations results 
showed that the interface between Hf02 and Si is found to 
be stable with respect to the formation of silicides whereas 
the Zr0 2 /Si interface is not[5] . The stability of the Hf0 2 /Si 
interface make the use of this oxide even more promising. 

The other important group IVA transition metals, such 
as Ti and Zr, which show very similar behavior as Hf, 
with oxygen are modeled recently by Waldner and Eriks- 
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son [6] and Wang et al. [7] respectively. All these systems 
commonly have a wide oxygen solubility ranges in the hep 
phase, up to 33 at.%(Ti-0), 29 at.%(Zr-0), and 20 at.%(Hf- 
O) at room temperature, as derived from higher tempera- 
ture measurements. 

In the present work, the Hf-O system has been modeled 
with the existing experimental data and first-principles 
calculations results. Afterwards, combining with the ther- 
modynamic parameters of the Hf-Si and the Si-0 binary 
systems, the thermodynamic description of the Hf-Si-0 
ternary system is obtained, and the stability diagrams 
pertinent to thin film processing, such as the Hf02-Si02 
pseudo-binary the isopleth of Hf02-Si, and isothermal 
sections are calculated. 

2. Experimental data 

2.1. Phase diagram data 
2.1.1. Hf-O 

Many investigations have been conducted to clarify the 
phase diagram of the Hf-O system[8, 9, 10, 11]. The main 
issues regarding the phase diagram of the Hf-O system in- 
clude: the extent of the a-Hf solid solution, the congru- 
ent melting of Hf02 , and the allotropic transformations of 
Hf0 2 . 

The phase diagrams suggested by Rudy and Stecher [8] 
and Domagala and Ruh [9] arc quite similar to each other, 
except the formation of the /3-Hf phase and the eutectic re- 
action around 37 at.% oxygen of liquid. For the hafnium 
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rich-side, Rudy and Stecher [8] proposed a eutectic reaction 
at 2273K, while Domagala and Ruh [9] suggested a peri- 
tectic reaction around 2523K. Similar to other group IVA 
transition metals, such as Ti, and Zr, it is strongly believed 
that Hf also has a peritectic reaction at the Hf-rich side. 
Another minor disagreement between these two phase dia- 
gram determinations is the eutectic reaction, Liquid —> a + 
H1O2. Rudy and Stecher suggested composition of oxygen 
in the liquid as 40 at. % O and 2453±40K, while Domgala 
and Ruh proposed 37 at.% O and -2473K. 

For the a-Hf solid solution, these two works[8, 9] are in 
quite good agreement with each other. Rudy and Stecher 
found that a-Hf dissolves up to 20.5 at.% oxygen at 1623K 
and that the solubility range is almost independent of tem- 
perature; this shows consistency with observations by Do- 
magala and Ruh and they quoted a solubility of oxygen in 
a-Hf of 18.6 at.% at 1273K. 

Ruh and Patel [11] proposed a tentative phase dia- 
gram for the Hf02-rich portion of the Hf-Hf02 system 
on the basis of metallographic data. They suggested the 
existence of solid solution regions for both cubic and 
tetragonal phases deviated from the stoichiometric com- 
position of HfC>2. Since the important phases in the Hf-0 
system are the polymorphs of the HfC>2 phases, i.e. mon- 
oclinic, tetragonal, and cubic, they have been studied 
extensively [12, 13, 14, 15, 16]. Allotropic transformations 
of the HfC>2 phase have been well summarized by Wang 
et al. [17]. 

The suggested phase diagram of the Hf-0 system by Mas- 
salski [18] is shown in Fig. 1 . 
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Fig. 1. Proposed phase diagram of the Hf-O system from Massalski 
[18]. 



2.1.2. Hf-Si-0 

Not many studies have been conducted regarding the 
phase stabilities of the Hf-Si-0 ternary system. Speer and 
Cooper [19] reported a ternary compound, Hafnon, with 
the chemical formula of HfSi04 and the crystal structure 
as IAi/amd. One of the most important phase diagrams of 



the Hf-Si-0 system is the Hf02-Si02 pseudo-binary that 
includes HfSi04 since the phase stabilities of this pseudo- 
binary are pertinent to the processing of the dielectric thin 
film. Parfenenkov et al. [20] determined the melting of 
HfSi0 4 at 2023±15K. 

2.2. Thermochemical data 

As discussed in the introduction, the Hf-0 system 
has a wide range of oxygen solubility in the hep phase. 
Hirabayashi et al. studied order/disorder transformation of 
interstitial oxygen in hafnium[21] around 10 ~ 20at.% by 
electron microscopy and neutron and X-ray diffractions. 
This work revealed that two types of interstitial superstruc- 
tures are formed in the hypo- and hyper-stoichiometric 
compositions near HfOi/6 below 700K which have i?3 and 
P31c symmetries, respectively. 

For the completely disordered hep phase at high tem- 
perature, Boureau and Gerdanian [22] measured the par- 
tial molar enthalpy of solution of oxygen in a-Hf solid so- 
lution at 1323K as a function of oxygen content using a 
Tian-Calvet-type microcalorimctcr. Previously, it was al- 
most impossible to measure the extremely low oxygen pres- 
sure in equilibrium with hafnium-oxygen solutions. There- 
fore derivation from the second law of thermodynamics 
was the only way to acquire thermodynamic information 
of solid solution phases[23, 24]. The major difficulty of di- 
rect measurement is making sure that all the hafnium sur- 
face is accessible at the same time to oxygen. Boureau and 
Gerdanian improved the accuracy of measurement by solv- 
ing the geometrical effect of specimen and oxygen contact. 
The observed phase boundary of a-Hf in this work is con- 
sistent with that of previous phase diagram studies [8, 9] as 
O/Hf=0.255. 

3. First-principles calculations 

3.1. Methodology 

The Vienna Ab initio Simulation Package (VASP)[25] 
was used to perform the Density Functional Theory (DFT) 
electronic structure calculations. The projector augmented 
wave (PAW) method [26] was chosen and the general gra- 
dient approximation (GGA)[27] was used to take into ac- 
count exchange and correlation contributions to the hamil- 
tonian of the ion-electron system. An energy cutoff of 500 
eV was used to calculate the electronic structures of all 
the compounds. 5,000 fc-points per reciprocal atom based 
on the Monkhorst-Pack scheme for the Brillouin-zone sam- 
pling was used. The fc-point meshes were centered at the T 
point for the hep calculations. 

3.2. Ordered phases 

The ordered structures of the Hf-Si-0 system calculated 
in this work can be categorized into three groups. First, 
pure elements, i.e. hep hafnium and diamond silicon are 
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calculated for the reference states. Second, hypothetical 
compounds, the end- members of the a, (3 solid solutions 
(HfOo.5 and HfC^), are also calculated. The stable com- 
pounds, monoclinic HTO2, quartz SiC>2, and the ternary 
compounds, HfSiC>4 are calculated as well. The calculated 
results of the ordered structures are listed in Table 1 . The 
enthalpy of those compounds are calculated from Eqn. 1: 

AH Hf. vStu O z =H(HfxSty0z) _ 

- -i^H(Si) - -J^-HiO) (1) 

x+y+z \ / x+yA-z V 1 \ / 

where H corresponds to the enthalpies of the compound 
and reference structures. The reference states for Hf and Si 
were the hep and diamond structures, respectively. In the 
case of condensed phases, the effects of lattice vibrations 
and other degrees of freedom (i.e. electronic, magnetic) can 
be neglected at low temperatures. Moreover, due to their 
rather small molar volumes, the PV contributions to their 
enthalpies (H = U + PV) can also be neglected. Thus 
their enthalpies H, can be replaced by the calculated first- 
principles total energies at OK in Eqn. 1. For the oxygen gas, 
the selected reference state was diatomic oxygen, O2. In this 
case, the contributions due to vibrational and translational 
degrees of freedom, as well as the PV work term, the molar 
volume of O2 is much larger than that of the condensed 
phases and cannot be neglected. In the following section, 
the determination of the correct reference state for O2 will 
be briefly discussed. 

3.3. Oxygen calculation 

The selected reference state for oxygen corresponds to 
the diatomic molecule at 298K. Two oxygen atoms were 
placed in a 'big box' to model the oxygen molecule gas and 
completely relaxed to find the lowest energy configuration. 
In order to properly account for the net magnetic moment 
of this molecule, spin polarization was considered. It was 
necessary to take into account the contributions of vibra- 
tional, rotational and translational degrees of freedom at 
finite temperature. Under the harmonic oscillator-rigid ro- 
tor approximation at temperatures greater than the char- 
acteristic rotational temperature — 2.07K for O2 — the in- 
ternal energy of the O2 molecule is given by [28] : 

where ©„ is the characteristic vibrational temperature. The 
first term in Eqn. 2 corresponds to the contributions due 
to translational and rotational degrees of freedom and the 
second and third terms correspond to vibrational contribu- 
tions. The last term, Eo, corresponds to the energy of the 
ground electronic state at OK. By including the PV = ksT 
term of an ideal, non-interacting gas in Eqn. 2, the enthalpy 
of the diatomic O2 molecule can be obtained. The differ- 
ence between the ground state electronic energy and the 
enthalpy of O2, per atom, is given by: 



H(T)0>-EO> = k jf Q + ^ + (3 ) 

In the case of O2 , the characteristic vibrational tempera- 
ture, 6„ is 2256K [28]. At 298K, the value of H(T)° 2 -E° 2 
is +0.0936 eV/atom or +9.03 kJ/mol-atom. 

3. 4- Disordered phases: Special quasirandom structures 

In order to calculate the enthalpies of mixing of oxy- 
gen and vacancy for both hep and bec phases in the Hf-0 
system, special quasirandom structures (SQS)[29] arc used 
which are supercells with correlation functions as close as 
possible to those of a completely random solution phase. 

In order to introduce the special quasirandom structure, 
it is convenient to understand the concept of correlation 
functions. In any atomic arrangements, the geometrical cor- 
relation between atoms in the structure can be defined as 
to whether it is ordered or disordered. For a completely dis- 
ordered structure, the surrounding environment of an atom 
at any given sites should be the same as all the other lattice 
sites. Correlation functions, \\ k , are defined as the prod- 
ucts of site occupation numbers of different figures, k, such 
as point, pair, triplet (when k = 1,2,3,...) and so forth. 
The correlation functions of each figure can be grouped 
together based on the distance from a lattice site as m th 
nearest-neighbors. 

For a binary system, a spin variable, a = ±1, can be as- 
signed to different types of atomic occupations and their 
products represent the correlation functions of binary al- 
loys. The correlation function of a random alloy is simply 
described as (2x — l) k in the A\- X B X substitutional binary 
alloy, where x is the composition. Once such a superccll sat- 
isfies the correlation function of a target structure, it can 
be easily transferred to other systems by simply switching 
the types of atoms in the structure. 

The major drawback of SQS is that the concentration 
which can be calculated is typically limited to 25, 50, and 
75 at.% since the correlation functions for completely ran- 
dom structures other than those three compositions are al- 
most impossible to satisfy with a small number of atoms. 
In principle, one can find a bigger supercell which has bet- 
ter correlation functions than smaller ones; however, such 
a calculation requires expensive computing. On the other 
hand, three data points from SQS calculations can clearly 
indicate the mixing behavior of solution phases. Another 
disadvantage of SQS is that it cannot consider the long 
range interaction since the size of the structure itself is lim- 
ited. It is reported that SQS works well with a system where 
short range interactions are dominant [30, 31]. 

Three different compositions, i.e. yo = 0.25, 0.5, and 0.75 
with yo representing the mole fraction of oxygen in the 
hep and bec interstitial sites, were considered and only two 
structures were generated in both phases since the struc- 
tures of yo = 0.25 and 0.75 are switchable to each other. 
For the a solid solution, the total number of lattice sites con- 
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Tabic 1 

First-principles calculation results of pure elements, hypothetical compounds (a, /3-Hf), and stable compounds (HfC>2, S1O2, and HfSiC>4). 
By definition, AH? of pure elements are zero. Reference states for all the compounds arc SER. 

Phases Space Lattice parameters (A) Total energy AH? 

Group a b c (eV/atom) (kj/mol-atom) 

HCP.A3 (Hf) P6 3 /rnmc 3.198 3.198 5.053 -9.8320 

Diamond_A4 (Si) Fd3m 5.468 5.468 5.468 -5.4315 

Gas (O2) -4.7936 

a-Hf(HfO . 5 ) P3ml 3.225 3.225 5.150 -9.9718 -175.511 

/3-Hf(Hf0 3 ) ImZm 4.364 4.364 4.364 -7.7253 -161.308 

Monoclinic(Hf0 2 ) a P2i/c 5.135 5.194 5.314 -10.2101 -360.563 

Quartz (Si0 2 ) P3 2 21 5.007 5.007 5.496 -7.9581 -284.809 

HfSiO-4 IAx/amd 6.616 6.616 6.004 -9.1024 -324.453 

-1.769 b 



a /3=99.56° 

b Reference states are monoclinic (Hf02) and quartz (Si02). 

sidered were 24, 36, and 48. Since only the mixing between 
oxygen and vacancies are considered to generate a SQS for 
the hep and bec phases, the hafnium ions are excluded from 
the correlation function calculations. Therefore, the total 
number of oxygen and vacancies are 8, 12, and 16, respec- 
tively. For the f3 solid solution, the total number of lattice 
sites that were considered for the mixing of oxygen and 
vacancy were 12 and 24 with total number of sites being 
16 and 32, respectively. The complete descriptions of the 
SQS's for a and (5 solid solutions are listed in Tables 2 and 
3, and their correlation functions are given in Tables 4 and 
5. Finding bigger cells than these were prohibited by the 
limited computing resources. 

Table 4 

Pair and multi-site correlation functions of SQS-./V structures for a 
solid solution when the c/a ratio is ideal. The number in the square 
bracket next to m is the number of equivalent figures at the same 
distance in the structure. 



Oxygen % 






50 




75 


SQS-JV 


Random 


21 


36 


18 


Random 


48 


n 2 ,i[3] 














0.25 


0.25 


n 2 , 2 [i] 














0.25 


0.25 


H 2 ,3[3] 














0.25 


0.25 


n 2)4 [6] 








-0.16667 





0.25 


0.20833 


n 2 , 6 [3] 








-0.11111 





0.25 


0.25 


n 2 ,6[3] 





-0.16667 0.11111 





0.25 


0.41667 


n 3 ,i[2] 








-0.33333 


-0.25 


0.125 


0.25 


n 3 ,i[6] 








0.11111 


-0.08333 


0.125 


0.08333 


n 3 , 2 [2] 








0.33333 


0.25 


0.125 


0.25 



The generated SQS's were either fully relaxed, or relaxed 
without allowing local ion relaxations, i.e. only volume for 
bec and volume as well as c/a ratio for hep were optimized. 
Theoretically, all the first-principles calculations should be 
fully relaxed to find the lowest energy configurations. How- 
ever, the structure should lie on the energy curve vs. geo- 



Table 5 

Pair and multi-site correlation functions of SQS-N structures for (3 
solid solution. The number in the square bracket next to il^. m is the 
number of equivalent figures at the same distance in the structure. 



Oxygen 


% 


50 




75 


SQS-7V 


Random 


16 


32 


Random 


32 


22,1 [6] 











0.25 


0.25 


5 2 , 2 [12] 











0.25 


0.25 


n 2 , 3 [i2] 











0.25 


0.25 


n 2 , 4 [6] 





0.16667 





0.25 


0.33333 


n 2 , 4 [3] 











0.25 


0.33333 


n 2 , 5 [24] 





-0.29167 





0.25 


0.25 


n 2 , 6 [24] 





-0.08333 





0.25 


0.25 


n 3 ,i[i2] 





-0.16667 





0.125 


0.16667 


n 3 , 2 [8] 











0.125 





n 3 , 3 [48] 





0.08333 0.08333 


0.125 


0.125 



metrical degree of freedom of the same phase. If the fully 
relaxed final structure does not have the same crystal struc- 
ture as initial input, it is not the phase of interest any 
longer. Thus it is necessary to force the structure to keep 
its parent symmetry. 

The calculated results of a and (3 solid solution phases 
are listed in Table 6 and 7. More detailed discussion of 
calculation results of special quasirandom structures and 
optimized results are found in a later section. 

4. Thermodynamic modeling 

11. Hf-0 

Seven phases are modeled in the Hf-0 system: hep, bec, 
ionic liquid, gas, and three polymorphs of HIO2: mono- 
clinic, tetragonal, and cubic phases. Detailed discussions of 
individual phases arc given below. 



Dongwon Shin et al. / Computer Coupling of Phase Diagrams and Thermochemistry (2008) 1-0 



5 



Table 2 

Structural descriptions of the SQS--/V structures for the a solid solution. Lattice vectors and atom/vacancy positions are given in fractional 
coordinates of the supercell. Atomic positions are given for the ideal, unrelaxed hep sites. Translated Hf positions are not listed. Original Hf 
positions in the primitive cell are (0 0) and ^ ^). 



Oxygen % 
SQS-iV 



Lattice vectors 



O 



Va 



V-3 
1 

3 

1 3 
_ 2 

3 



21 






36 







f 1 


-1 


-2 




-2 


-2 









-3 


i 


3 


2 


-2 


3 


4 


3 




-1| 


3 
4 


X 3 


-2 


"If 


3 
4 


-22 

^3 


— 1 


"If 


-If 


2 
3 


-2 






x 3 


-4 






2 
3 


-4 



18 



75 
48 







\1 



1 

3 " 

3 L ' 
1 

3 " 

i 

3 " 
1 

3 " 



^ -It 



3± 



2 A 

-3 

3 3 



-1J 




1 

3 


-if 


*3 


-3i 

J 3 


1 

4 


1 


3 


2 




1 


3 


4 


3 


4 


1 

3 


3 
4 


^3 


"4 






3 


2 


-1| 


i 


4 


3 


4 






X 3 


-3i 
^3 








2 
3 


-3i 

J 3 





2 t 

l 

3 



1 


3 


1 


3 


4 


3 


1 


3 


1 


3 


4 


3 


"IS" 


3 
4 


13- 


2 


3 




3 


4 


2 


3 




3 


4 




1 

3 


-If 




1 

3 


-If 




"If 


"If 





Table 6 

First-principles calculations results of a-Hf special quasirandom 
structures. FR and SP represent 'Fully Relaxed' and 'Symmetry 
Preserved', respectively. Oxygen atoms are excluded for the symme- 
try check. 



Oxygen Atoms Symmetry Space Total Energy 
% Hf O Va Group (eV/atom) 



(kj/mol-atom) 



25 


32 


1 


12 


FR 


P63 /mmc 


-9.9140 


-61.9280 




32 


1 


12 


SP 


P&z/mmc 


-9.9072 


-61.2715 


50 


16 


1 


1 


FR 


PQ^/mmc 


-9.9590 


-109.481 




16 


4 


i 


SP 


P63 /mmc 


-9.9444 


-108.075 




21 


6 


6 


FR 


P63 /mmc 


-9.9560 


-109.187 




21 


6 


6 


SP 


PGs/mmc 


-9.9418 


-107.818 




32 


8 


8 


FR 


P63 /mmc 


-9.9564 


-109.230 




32 


8 


8 


SP 


P63 /mmc 


-9.9413 


-107.768 


75 


32 12 


4 


FR 


PG^/mmc 


-9.9755 


-146.428 




32 12 


4 


SP 


PG^/mmc 


-9.9541 


-144.356 



Tabic 7 

First-principles calculations results of /3-Hf special quasirandom 
structures. FR and SP represent 'Fully Relaxed' and 'Symmetry 
Preserved', respectively. Oxygen atoms are excluded for the symme- 
try check. 



Oxygen 


Atoms 


Symmetry Space Total Enerj 


;y AH mix 


% 


Hf O Va 




Group 


(eV/atom 


(kj/mol-atom) 


25 


8 6 18 


FR 


Pm 


-9.9280 


-227.299 




8 6 18 


SP 


Pm3m 


-9.0165 


-139.345 


50 


4 6 6 


FR 


C2/m 


-10.0088 


-315.520 




4 6 6 


SP 


Im3m 


-8.5718 


-176.866 




8 12 12 


FR 


Pm 


-9.9636 


-311.160 




8 12 12 


SP 


Im3m 


-8.6081 


-180.370 


75 


8 18 6 


FR 


Pm 


-9.4727 


-307.100 




8 18 6 


SP 


ImUm 


-8.1752 


-181.908 



4.1.1. HCPandBCC 

It has been reported that the stable solid phases of group 
IVA transition metals (Ti, Zr, and Hf), the hep and bec 
phases, dissolve oxygen interstitially into their octahedral 
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Table 3 

Structural descriptions of the SQS-./V structures for the (5 solid solution. Lattice vectors and atom/vacancy positions are given in fractional 
coordinates of the supercell. Atomic positions are given for the ideal, unrelaxed bec sites. Translated Hf positions are not listed. The original 
Hf position in the primitive cell is (0 0). 



Oxygen % 
SQS-iV 



50 



16 



32 



75 
32 





f 


0.5 


0.5 


1.5 > 




(- 


-1 





Lattice vectors 




-0.5 


1.5 


0.5 


| 
I 







1 




{- 


-0.5 - 


-1.5 - 


-0.5 J 




\ 





-2 













0.5 






1 







-0.5 


1 


1 






1 


0.5 




-0.5 





1 




-0 


5 







o 


—0.5 


1 






1 


0.5 







1 


1.5 




-0 


5 


0.5 







0.5 


1.5 






1 


-1 














1 - 


-0.5 














1 - 


-0.5 












-0 


5 

















1 - 


-0.5 












-0 


5 - 


-0.5 












-0 


5 


-1 



Va -0.5 


1 0.5 


-0.5 


-1 


-2.5 


-1 2 1.5 


-0.5 


-0.5 


-1 


-1 


-2.5 


-1 1.5 1.5 


-0.5 


-1 


-1 


-1.5 


-2.5 


-0.5 1 0.5 


-0.5 


0.5 


-0.5 





-0.5 


-1 1.5 2 





0.5 1 


-1 


-1.5 


-2 


-0.5 0.5 2 


-0.5 


0.5 1 


-0.5 


-1.5 


-2 


-0.5 2 2 






-0.5 


-1 


-3 








-1 





-1.5 








-1 


-0.5 


-1.5 








-0.5 





-1 








-0.5 


-1 


-1.5 








-0.5 


-0.5 


-1 





sites [32]. The solid solutions of Hf-0 are modeled with the 
two sublattice model, with one sublatticc occupied only by 
hafnium and the other one occupied by both oxygen and 
vacancies: 



where c corresponds to the ratio of interstitial sites to nor- 
mal sites in each structure. For the hep phase, the ratio 
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is derived as c = | from the pure hep sublattice model, 
consistent with the previous thermodynamic modeling of 
Ti-0[6] and Zr-0[33]. For the bcc phase, the stoichiometric 
ratio c is equal to 3 as in the conventional bcc phase. 

The Gibbs energies of the hep and bcc phase can be 
described as: 



n HCP.BCc „, on _i_ „, On 

+ cRT( yo ln(yo) + Wa HWa)) + G% 

(4) 

k 

=yoyva^2 k L (Hf: oy a) {T){yo - Wa) k (5) 

fc=0 

where Gnfo c is the standard Gibbs energy of the hypo- 
thetical oxide HfO c , which is one of the end members that 
establishes the reference surface for this model; °GHfVa c 
corresponds to the standard Gibbs energy of the pure 
bcc and hep phases and the chemical interaction for oxy- 
gen and vacancies in the second sublattice is given by 
yoyvaJ2 kL (Hf-.o,Va)(T)(yo ~ yva) k - This is identical to 
a Redlich-Kistcr polynomial [34]. 

4-1-2. Ionic liquid 

The liquid phase region goes from pure liquid hafnium to 
stoichiometric liquid Hf02 . The ionic two-sublattice model 
is used for the liquid phase [35]: 

(C+ Vi )p(Ap,Va,B° k ) Q 

where Cf Vi corresponds to the cation, i, with valence, +i>;; 
Aj to the anion, j, with valence, —vy, Va arc hypotheti- 
cal vacancies added for electro-neutrality when the liquid 
is away from stoichiometry, having a valence equal to the 
average charge of the cation, Q; and B\ represents any 
neutral component dissolved in the liquid. The numbers of 
sites in the sublattices, P and Q, are varied in such a way 
that electro-neutrality for all compositions is ensured with 
y being the site fractions: 

3 

Q = ^v l y Ci (7) 

i 

For the particular case of the Hf-0 system, this two- 
sublattice model can be further simplified: 

(Hf+%_ 2yo _ 2 (0-\Va-% 
The Gibbs energy expression for this system is: 



7 

G Ionic liquid „, 0/~iL 
ra -yPf+iyO- 2 (jr (J?/+ 4 ) 2 (ya-*) 4 

+ 4j/_H-/+ 4 yya- 4 °G(_H-/+4)(y a -4) 

+ RT(A - 2y -2)(y HS+i Hy Hf+4 )) 

+ -R r ( 4 )(j/o- 2 Myo-O + yva-i ln(yy a - 4 )) 
+ GZ (8) 

k 

=te/+ 4 2/o- 2 yya- 4 ^ k L(Hf+ i :0- 2 ,va-^)(T) 

k=0 

X (Z/O- 2 ~yVa-*) k (9) 

where G^ f y. +4 ^^ a _4^ corresponds to the standard Gibbs 

energy for two moles of liquid -ff/02; GtHf+^iVa- 4 -) lS 
the standard Gibbs energy for pure hafnium liquid and 
fe L(fly+4.Q-2.y a -4) corresponds to the excess chemical in- 
teraction parameters between oxygen and vacancies in the 
second sublattice. 

4.1.3. Gas 

To describe the oxygen-rich side of the Hf-0 system, the 
gas phase was included in the calculation. The ideal gas 
model was used and the following six species were consid- 
ered: 

(0,0 2 ,0 3 ,Hf, HfO,Hf0 2 ) 
The Gibbs energy of the gas phase can be described as: 

G,?r = yo( G o + RT\n(P)) + yo 2 ( G O2 + BTIn(P)) 
+ yo 3 (°Go 3 +RT\n(P)) 
+ yHf( a G Hf +RT\n(P)) 
+ yHfo{°G H fo + RTHP)) 
+ y Hf o 2 ( G Hf o 2 +RTln(P)) 
+ RT(yo ln(yo) + yo 2 Mvo 2 ) + Vo 3 ln(yo 3 )) 
+ RT{ynf ln(2/H/) + VHfo ^(y Hf0 ) 
+ y H fo 2 MvHfo 2 )) (10) 

with y being the mole fraction of species in the gas phase. 

4-1-4- Polymorphs of B.jO% 

Thermodynamic descriptions of three polymorphs of 
Hf02 have been obtained from the SSUB database[36]. For 
simplicity, all three phases are modeled as line compounds 
and the transformation temperatures for monoclinic — » 
tetragonal — > cubic — > liquid are 2100, 2793, and 3073K, 
respectively. 

4.2. Si-0 

The Si-0 system has been modeled by Hallstedt [37] with 
an ionic liquid model. Three different polymorphs of silicon 
dioxides, quartz, tridymite, and cristobalite are included in 
the system. The Si-0 phase diagram is given in Fig. 2. 

4.3. Hf-Si 

The Hf-Si system has been extensively studied and mod- 
eled by Zhao et aZ[38]. Six intcrmetallic compounds, IH^Si, 
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Fig. 2. Calculated Si-O phase diagram from Hallstedt[37]. 

Hf 5 Si 3 , Hf 3 Si 2 , Hf 5 Si 4 , HfSi, and HfSi 2 are present. The 
calculated phase diagram of Hf-Si system is given in Fig. 3. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
Mole Fraction, Si 

Fig. 3. Calculated Hf-Si phase diagram from Zhao et ai[38]. 



44. Hf-Si- o 

In order to be combined with the Hf-0 and Si-0 systems, 
the liquid phase of the Hf-Si system was converted to an 
ionic liquid in the present work. Hf +4 and Si +4 are in the 
first sublattice of the ionic liquid phase and vacancies have 
been introduced into the second sublattice for the electro- 
neutrality The interaction parameters for the liquid phase 
from Zhao et al. [38] were used for the mixing of Hf +4 and 
Si +4 in the first ionic liquid sublattice. 



The ternary compound HfSi04 has been introduced. Due 
to the lack of experimental data, HfSi04 has been modeled 
as a stoichiometric compound. 

5. Results and discussion 

Based on the existing experimental data and the first- 
principles calculations results, the model parameters for the 
Hf-0 and Hf-Si-0 systems are evaluated. The PARROT 
module in the Thermo-Calc software has been used [39]. 

The enthalpy of mixing of the hep and bec phases calcu- 
lated from the present thermodynamic modeling arc shown 
together with the result of first-principles calculations in 
Fig. 4. For the a solid solution, they agree well with each 
other. As shown in Table 6, all first-principles calculation 
results retained their original symmetry as hep and the cal- 
culation results are quite well converged with respect to the 
size of SQS. On the other hand, results for the (3 phase are 
not as good as those for the a phase. The difference be- 
tween fully relaxed and symmetry preserved calculations is 
140kJ/mol-atom at most and all the fully relaxed calcula- 
tions could not maintain the bec symmetry. 
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Fig. 4. First-principles calculations results of hypothetical com- 
pounds (HfOo.5 and HfOa) and special quasirandom structures for a 
and j3 solid solutions with the evaluated results. Reference states for 
Hf of a and /3 solid solutions are given as hep. Fully relaxed calcu- 
lations of fi solid solution have been excluded from this comparison 
since the calculation results completely lost their bec symmetry. 

Such first-principles calculations results can be validated 
by comparing their calculated crystallographic information 
with experimental measurements. In Fig. 5, the calculated 
lattice parameters of a-Hf, both a and c, are compared with 
experimental data and quite satisfactorily agree with each 
other. 

In Fig. 6, the Hf-rich side of the Hf-0 phase diagram 
is shown. The congruent melting of a-Hf and the peritec- 
tic reaction arc reproduced correctly The calculated phase 
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Fig. 5. Calculated lattice parameters of a-Hf with experimental 
data[8, 9, 24, 43]. Scale for a-axis is left and for c is right. 

region shows quite good agreement with the X-ray phase 
identification results of Domagala and Ruh [9] . 
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Fig. 6. Calculated Hf-rich side of the Hf-O phase diagram with 
experimental data from Domagala and Ruh [9], 

The calculated partial enthalpy of mixing of oxygen in 
the a-Hf is shown in Fig. 7 and compared with experimen- 
tal data[22] . As mentioned before, the accuracy of measure- 
ment has been improved compared to the previous experi- 
ments. However, it is still quite difficult to measure the low 
oxygen pressure. 

The calculated phase diagram of the entire Hf-0 system 
is shown in Fig. 8 with the gas phase included. 

In the present work, the ternary liquid phase is extrap- 
olated from the binaries. The enthalpy of formation of 
HfSiC>4 is obtained from first-principles calculations and 




0.05 0A0 0A5 020 0.25 0.30 
Mole Fraction, O 



Fig. 7. Calculated partial enthalpy of mixing of oxygen in the ct-Hf 
with experimental data[22] at 1323K. 
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Fig. 8. Calculated Hf-O phase diagram. 



the entropy of formation is evaluated from its preitectic re- 
action, Liquid + Hf0 2 (Monoclinic) -> HfSi0 4 , at 2023K as 
reported by Parfcncnko et al [20] . The pseudo-binary phase 
diagram of Hf02-Si02 is calculated and shown in Fig. 9. 

As discussed in the introduction, Hf02 is a promising 
candidate to replace SiC>2 as the gate dielectric in CMOS 
transistors due to its high dielectric constant and compati- 
bility with Si in comparison with Zr02_5]. In general, dur- 
ing the fabrication of such devices, the films are subjected 
to temperatures around 1273K for a short period of time [4] . 
Thus, it is quite essential to understand the thermodynamic 
stability of Hf0 2 /Si0 2 /Si. 

From the evaluated thermodynamic database of the 
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Fig. 9. Calculated HfC>2-SiC>2 pseudo-binary phase diagram. 

Hf-Si-O, the isothermal sections of the Hf-Si-0 system 
can be readily calculated to study the stability of the 
HfC>2/Si interface. Two different temperatures are selected 
for the calculations and they arc 500K for low tempera- 
ture processing, such as mist deposition method and rapid 
thermal processing [40] and 1000K, typical temperature 
for the epitaxial growth of oxides deposition. Calculated 
isothermal sections of the Hf-Si-0 system at 500K and 
1000K are shown in Fig. 10(a) and 10(b), respectively. 
The two three-phase regions, HfSi04+Hf02+Hf2Si and 
HfSi04+diamond+Hf2Si, in the 500K isothermal section 
should be noticed with respect to the stability of HfC^/Si 
interface. Since those regions are intersected by the line 
connecting Hf02 and Si, HfSi2 can be found in the fabrica- 
tion of polySi/Hf02 gate stack Metal Oxide Semiconductor 
Field Effect Transistor (MOSFET) on bulk Si at 500K, 
while the 1000K calculation result clearly shows that Hf02 
is stable with the Si substrate. 

A vertical cross section of isothermal sections, the iso- 
pleth between Hf02-Si is calculated in order to investigate 
the stability range of HfSi2 in the HfO^/Si interface and 
is given in Fig. 11. It should be noted that HfSiCU at the 
low temperature range is zero amount. The calculated re- 
sult shows that HfSi2 becomes stable below 543. 53K. This 
result is in agreement with the experimental observation 
from Gutowski et al[5]. Their HTO2 film was deposited at 
823K and then annealed at 1023K without the formation 
of any silicides. 

It should be emphasized that the thermodynamic stabil- 
ity of HfSi2 in the HfC^/Si interface depends on the for- 
mation energy of HfSi04 based on the isothermal sections 
and the isopleth calculations. The enthalpy of formation 
for HfSi04 is calculated from first-principles calculations 
since there is no experimental measurement. To further il- 
lustrate this, the reference states of the enthalpy of forma- 



5000 
4500 
4000 
3500 
u - 3000 
I 2500 

a 2000 

H 

1500 

1000 H 

500 





L1+L2 
^HfCtyt) 

Gas+Li+Hf0 2 (t) Li+L 2 +Hf02(m) Li+L: 



=1= 



Li+Hf0 2 (m)+HfSi04 



Li+HfSiO. 



1 



Hf02(m)+diamond[+L:] 
Hf02(m)+diamond[+HfSi04] 



L1+L2 
+HfSi04 

543.53 



HfSi04+Hf02(m) I HfSi04+diamond 
+HfSi2 +HfSi2 
I I 1 I 1 I I I I 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

Mole Fraction, Si 



Fig. 11. Calculated isopleth of Hf02-Si at 1 atm. Hafnium dioxide is 
left and silicon is right. Polymorphs of Hf02, monoclinic, tetragonal, 
and cubic, are given in parentheses. The phases in the bracket are 
zero amount. 

tion for HfSi04 arc defined as the two binary metal oxides 
(See Eqn. 11 and Table 1). 



AHf fSi0i = H{HfSi0 4 ) - ±H{Hf0 2 ) - \H(Si0 2 ) 

(11) 

The current result of the HfSi04 calculation predicts that 
HfSi2 is stable up to 543. 53K. However, the uncertainty of 
the formation enthalpy of HfSi04, which originates from 
the density functional theory itself, is about ±1 kj/mol- 
atom[41]. Thus, the associated decomposition temperature 
of HfSi2 in the HfO^/Si interface varies from 381. 92K to 
669. 76K within the calculated uncertainty of HfSiC>4. The 
recent work from Miyata et al. [42] found the formation of 
nanometer-scale HfSi2 dots on the newly opened void sur- 
face produced by the decomposition of Hf02/Si02 films 
at the oxide/void boundary in vacuum. However, their re- 
sult cannot be directly compared with the current thermo- 
dynamic calculations due to the unknown oxygen partial 
pressure. 

All parameters for the Hf-Si-0 system are listed in Table 



6. Summary and conclusion 

The complete thermodynamic description of the Hf-Si- 
O ternary system is developed via the hybrid approach 
of first-principles calculations and CALPHAD modeling in 
the present work. In the Hf-0 system, special quasirandom 
structures have been generated to calculate the enthalpies 
of mixing of oxygen and vacancies in the a and (3 solid 
solutions. Calculated enthalpies of mixing of a-Hf are al- 
most identical to the model-calculated value whereas those 
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Fig. 10. Calculated isothermal section of Hf-Si-O at (a) 500K and (b) 1000K at 1 atm. Tie lines are drawn inside the two phase regions. The 
vertical cross section between Hf02 and Si is the isopleth in Figure 11. 



Table 8 

Thermodynamic parameters of the Hf-Si-O ternary system (in S.I. units). Gibbs energies for pure elements and gas phases are respectively 
from the SGTE pure elements database[44] and the SSUB database[36]. 

Phase Sublattice model Evaluated description 



Ionic 
liquid 



(Hf+^pCO-^Va^OG^l 1 ^ 



op monoclinic 
/(j Hf0 2 



252000 - 86.798T 



hep (Hf)i(0,Va)o. 5 



Q^Ionic liquid Q^Liquid 



Hf+ 4 :Va 
0y Ionic liquid 

Hf+-»:0- 2 ,Va : 
\ j Ionic liquid 

Hf+4:0- 2 ,Va : 
2y Ionic liquid 

Hf+4:0- 2 ,Va : 



Hi 

50821 



4.203T 
420485 - 133.300T 
30537 



Ophcp 
^HfiVa 

Ophcp _ 
^Hf:0 " 

Ot hc P 

Hf:0,Va 
l T hcp 



_ 0p hc P 
" ^Hf 
Op hep i 
°Hf H 



■0.5°Gg as 



271214 + 41. 560T 



-31345 
-6272 



bec (Hf)i(0,Va) 3 



HfSiQ 4 (Hf)i(Si)i(0) 4 



Opbcc 
^Hf:Va 

0("<bcc _ 
^Hf:0 " 

0t hc P 
^HfiO.Va 

0^.HfSiO 4 



_ Ofbcc 

Ophcp , oOpG 
^Hf +d 



: -981440 - 

(~* monoclinic 

Kj mo 2 



737857 - 
20.349T 

^ quartz 
°Si0 2 



268.540T 



10615 + 1.313T 



of /3-Hf show significant discrepancy. In the /3 phase, first- 
principles calculations could not retain its original symme- 
try as bec due to the strong interaction between the atoms 
in the structure. The calculated enthalpies of mixing from 
SQS's results are combined with the enthalpies of forma- 
tions of those hypothetical compounds calculated from the 
electronic structure calculations to derive the Gibbs energy 
of solid solutions in the Hf-0 system. 

In the total energy calculation of oxygen gas, vibrational, 
rotational and translational degrees of freedom are consid- 
ered. With the adjusted total energy of oxygen molecule, 
the enthalpies of formation from first-principles calcula- 
tions for both ordered and disordered phases showed good 
agreement with evaluated values. 



The Hf-0 system has been combined with the Hf-Si and 
the Si-0 systems to calculate the Hf-Si-O ternary system 
with the ternary compound HfSiO^ introduced from the 
first-principles calculations. From the Hf-Si-O thermody- 
namic database, phase stabilities pertinent to thin film pro- 
cessing such as Hf02-Si02 pseudo-binary, isothermal sec- 
tions, and isopleth have been calculated. The thermody- 
namic calculation results show that the Hf02/Si interface 
is stable above 543. 53K, which agrees with previous exper- 
imental results. However, due to the uncertainty of HfSi04 
formation energy from first-principles, the stability of HfSi2 
in the Hf02/Si interface is still in question and further ex- 
perimental investigation is required. 

It can be concluded that the thermodynamic properties 
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of solid phases can be obtained from first-principles calcu- 
lations not only for the ordered structures but also for the 
solution phases as long as one can find appropriate geo- 
metrical input for phases of interest. Special quasirandom 
structures for a substitutional solution phase is one exam- 
ple. However, one should notice that such a supercell only 
mimics short-ranged interaction as in metallic alloy sys- 
tems. As shown in the present work, SQS's can successfully 
describe the mixing behavior between oxygen and vacan- 
cies in the a solid solution where oxygen concentrations are 
relatively low, but for the oxygen-rich (5 phase such interac- 
tions between the electrons at the longer distance become 
important and lead to the collapse of its original structure 
as bec when the structure has been fully relaxed. 
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